Lyapunov Instability for a hard-disk fluid in equilibrium and nonequilibrium 
thermostated by deterministic scattering 
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We compute the full Lyapunov spectra for a hard-disk fluid under temperature gradient and 
shear. The system is thermalized by deterministic and time-reversible scattering at the boundary. 
This thermostating mechanism allows for energy fluctuations around a mean value which is reflected 
by only two vanishing Lyapunov exponents in equilibrium and nonequilibrium. The Lyapunov expo- 
nents are calculated with a recently developed formalism for systems with elastic hard collisions. In 
a nonequilibrium steady state the average phase-space volume is contracted onto a fractal attractor 
leading to a negative sum of Lyapunov exponents. Since the system is driven inhomogeneously we 
do not expect the conjugate pairing rule to hold which is confirmed numerically. 
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q'. I. INTRODUCTION 

u. Transport of energy and momentum is a central problem in nonequilibrium statistical mechanics, but so far most 
of our knowledge is confined to the macroscopic level. There is still a long way to go when it comes to understanding 
this phenomena in microscopic terms although significant progress has been made during the last years, with help 
from dynamical systems theory and computer simulations. In general, external forces are needed to drive a system out 
of equilibrium, but in order to prepare a nonequilibrium steady state the redundant energy has to be removed thus 
H ' preventing the system from heating up indefinitely. One way out is the introduction of thermostating mechanisms ■ 
| Both stochastic Q-[b§] and deterministic/time-reversible thermostats are in use, the latter having been introduced to 
t-H , remain close to Hamiltonian dynamics. 

Recently, an alternative thermostating mechanism [p~5| — [T^] acting via deterministic time-reversible boundary- 
scattering has been applied on a hard disk fluid to model heat and shear flow nonequilibrium steady states [fj"8| . 
The calculated transport coefficients have been found to be in agreement with the theoretical values obtained from 
kinetic theory, but only for special cases and in the thermodynamic limit the conjectured identity between exponential 
phase-space contraction and entropy production rate holds. In the present paper we investigate further the dynam- 
ical properties of this system by computing the full Lyapunov spectra and related quantities like the Kaplan- Yorke 
dimension or th e K olmogorov-Sinai entropy. In Sec. n] we briefly recapitulate the model and its thermostating mech- 
anism and Sec. Ill serves to outline the method used for computing the Lyapunov exponents [l9|p0[ |. The results are 
presented in Sec. IV and conclusions are drawn in Sec. 0. 



II. MODEL 

Consider a two-dimensional system of hard disks confined in a square box of length L with periodic boundary 
conditions along the x-axis, i.e., the left and right sides at x = ±L/2 are identified. The TV disks interact among 
themselves via elastic hard collisions, thus the bulk dynamics is purely conservative. In the following and in all the 
numerical computations we use reduced units by setting the particle mass m, the disk diameter a and the Boltzmann 
constant Ub equal to one. Now, denote with p l xl p l y and with pf, p? the tangential and normal momentum of a disk 
before and after a collision with the wall. Then the scattering prescription is given as |l8| 

( V f v f) = { T ~ lo M°T{P^P% Pl>0 ( ) 

{Px:Py) \ U oT - 1 CM' 1 oT (pl,pl), pl<0, {! 

where T : [0, oo) x [0, oo) — > [0, 1] x [0, 1] is the invertible map 

(CO - T (p x ,p v ) = (erf (|jfc|/V2?) , cxp (-p 2 y /2T)) (2) 
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and M. : [0, 1] x [0, 1] — > [0, 1] x [0, 1] is a two-dimensional, invertible, phase-space conserving chaotic map to be specified 
later. H(p x ,p y ) = (p x , —p v ) only serves to produce the right sign for the backward scattering. The parameter T plays 
the role of a temperature |15|-|lq] . Note that the colliding disk retains its tangential direction and that the scattering 
is reversible by construction. 

So far, the model has only been defined in equilibrium. In order to drive the system into a n oncqu ilibrium steady 
state (NSS) the collision rule only has to be modified appropriately, which will be done in Sec. IV B (see also [pi). 



III. METHOD 



Like all hard disk systems our model is chaotic in the sense that two nearby phase space trajectories diverge 
exponentially with time. This is mainly due to the dispersing action of the hard disk collisions in the bulk, but 
especially for small particle numbers we also expect a contribution to this divergence due the chaotic nature of our 
scattering mechanism. The average logarithmic divergence rate in phase space are described by the so-called Lyapunov 
exponents A;. Denote by T = {qi, q2, qw, Pi, P2, Pn} the 4N dimensional phase space vector for N disks. Then, 
the time evolution 

r(0) = $*[T(0)] (3) 

of an initial state T(0) consists of a smooth streaming which is interrupted by particle-particle and particle- wall 
collisions. Next, consider a satellite trajectory T s (t) initially displaced from the reference trajectory by an infinitesimal 
vector ^r(O). In a chaotic system | <5r (0) | is growing on average exponentially, thus rendering the system unpredictable 
for long times. Then there exists a complete set of linear-independent initial vectors {SFi(0) : I = 1, . . . ,4iV} and 
Lyapunov exponents defined as ]2l| ] 

A<=hmiln^% (4) 
t-oo t |«5T,(0)| V ' 

The A/, which we order according to Ai > A 2 > . . . > X^n, are independent of the coordinate system and the metric. 
The whole set of Lyapunov exponents is referred to as the Lyapunov spectrum. 

In Hamiltonian systems the Lyapunov exponents appear in pairs summing up to zero, Ai + XiN-i+i = for 
i = 1, ...,2N, due to the symplectic nature of the equations of motion. In a continuous dynamical system one 
Lyapunov exponent associated with the direction of the phase flow vanishes. Moreover, each conserved quantity leads 
to an additional vanishing Lyapunov exponent. The symmetry found in symplectic dynamical systems is lost when 
the system is driven to a nonequilibrium stationary state. However, for homogeneous driving the symmetry is replaced 
by the so-called conjugate pairing rule |22| saying that after excluding the vanishing exponents associated with the 
flow direction and the conservation of energy the remaining pairs, i.e. {Ai, \4n}, {A2, A4jv_i}, and so on, each sum 
up to the same negative value C. For inhomogeneously driven systems such as ours or the Chernov-Lebowitz shear 
flow model |ll],[l2|,^0|] , however, the symmetry is lost and no pairing rules exist. 

As can be seen from Eqs.(P^) the phase space volume is in general changed during each disk- wall collision. In 
equilibrium this averages up to zero, whereas in NSS the average phase space contraction rate is negative and is 
given by the sum of all Lyapunov exponents. Consequently, the phase volume shrinks continuously in NSS and and 
the phase-space distribution collapses onto a multifractal strange attractor. The fractal dimension of this strange 
attractor can be estimated with the conjecture of Kaplan- Yorke |23[ |, 

D KY = n+^± (5) 

where n is the largest integer for which Yl^—i A; > 0. Dky is the dimension of a phase space object which neither 
shrinks nor grows and for which the natural measure is conserved by the flow. 

For the calculation of the full Lyapunov spectrum we use a method worked out by Dellago et al. (l!]j2(| which is 
actually a generalization of the algorithm of Benettin et al. (24|] for smooth dynamical systems. The latter follows 
the time evolution of a reference trajectory and of a complement set of tangent vectors by solving the original and 
the linearized equations of motion, respectively. Periodic reorthonormalization prevents the tangent vectors from 
collapsing all into the direction of fastest growth. Averaging the logarithmic expansion and contraction rates of the 
tangent vectors then yields the Lyapunov exponents. For a hard disk system the free streaming is interrupted by 
impulsive collisions, either with another particle or the boundary. This certainly affects both the trajectory and the 
tangent space and has to be included in the calculation. The free streaming and the particle-particle collisions in the 
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bulk have been treated in Section III-D of reference [^9[ , where the same notation was used as in the present work. We 
refer to Eqs. (39)-(42) and Eqs. (68)-(73) of this article for explicit expressions of the particle-particle collision rules 
in phase space and tangent space, respectively. It remains to consider the particle-wall collisions and the follo wing 
lines are formulated in parallel to the respective treatment of Dcllago and Posch for the Chernov-Lebowitz model [E0| . 
In fact, the only difference lies in the 'scattering matrix' and its derivatives. 

If particle k collides with the walls its position remains unchanged whereas its momentum is changed according to 
the scattering rules Eq. (|l|). The collision map = M(r l ) in phase space becomes 

q\ q; for j = l,...,N (6) 

pj = p} for j^k (7) 

p{=C+(p 4 fc ), for pl>0 (8) 

p{=C-(p* fe ), for pl<0 (9) 

using the abbreviations C+ =2 1 o M. o T for the forward scattering and C~ = II o 3 1 o Ad^ 1 o T for the 
backward scattering (Eq. ([!])). 

In order to obtain the corresponding transformation for the tangent space vector ST at a particle-wall collision we 
assume that the collision takes place at phase point T at time r c . Then the satellite trajectory, displaced by the 
infinitesimal vector 5T, collides at a different phase point r + <5r c at a different time t c + 6t c . A linear approximation 
in phase space and time yields (l9) 



or 



dM 

— F(T)-F(M(r)) 



(10) 



where F is the right hand side of the equation of motion during the free streaming |l(| , and <9M / dT is the matrix of 
the derivatives of the full collision map with respect to the phase-space coordinates. Obviously, the delay time 6t c is 
a function of the phase point T l and of the tangent vector 5T l . For a disk- wall collision of the fcth particle the delay 
time St c is given by 

(Sq k ■ n) 

0T C =-- r. 11 

(Pfe/m • n) 

Here, n is the normal vector of the wall pointing into the simulation box. Since the scattering rules Eq. (|l]) for the 
momentum components is independent of the position of the particle, the matrix dlsA/dT has the form 

where 1 and are the 2N x 2N unit and zero matrices, respectively. dC ± (p l )/d(p z ) is the matrix of the derivatives 
of the outgoing momenta with respect to the incoming momenta and only the components of the colliding particle k 
are different from zero. From Eq. ([Io|) the following transformation rules for the tangent vectors can be deduced: 

Sqjj = 6q) for j^k (13) 
Spf = Sp) for j^k (14) 
Sq{ = 54 - (p[ - pl)Sr c (15) 

*d = ^S»l (16) 

Omitting for notational convenience the index k indicating the colliding particle we obtain from Eq. (Q) the following 
expressions for the 2x2 matrix (dC ± (Pk)/d(p l k )) a f3 — dp^/dp^, a,/3 G {x,y}: 

M = (DM) u exp [({p{f - {pl) 2 )/(2T)} , M = - (DM) 12 ^Lexp [((pif - ( P \) 2 )/(2T)\ (17) 

M = - (DM)ai ^ exp [((pff - (pl) 2 )/(2T)] , M = ( D - M ) 22 4 ex P Wv? ~ (^) 2 )/( 2T )] - ( 18 ) 
°Px -/>:, op y p j y 

for pi > 0. 
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Here, T)A4 denotes the matrix of the derivatives of the chaotic map M.. Eqs. ( |17|jl8[ ) are stated for positive tangential 
velocities, for negative tangential velocities M has only to be replaced by M. , see Eq. (0). 

Combining the free streaming with the transformation for the disk-disk and the disk-wall collisions, one is now able 
to follow the exact time evolution of the trajectory and of the tangent-space vector. 

IV. RESULTS 

Using the algorithm outlined in the previous section we are now able to calculate the full Lyapunov spectrum for 
our hard disk model with deterministic scattering at the boundary. As already mentioned we use reduced units by 
setting the particle mass m, the disk diameter a and the Boltzmann constant ks equal to unity. We define the 
number density by n = N/ L 2 . For simulation we use a collision-to-collision approach and neighbor lists |t]] . For an 
initial configuration the centers of the disks are positioned on a triangular lattice and the momenta are chosen from a 
Gaussian with zero mean. The total momentum is then set to zero and the momenta are rescaled to obtain the total 
kinetic energy Ekin = N(T U + T e i)/2, T u and Td being the imposed 'parametrical' temperatures of the upper and the 
lower wall, respectively. 

A. Equilibrium 

We now set both wall temperatures T Ul Td equal to one and compute the full Lyapunov spectra for a four-particle 
system at number density n — 0.2 using three different chaotic maps: 

M B {C0 = <Xj£/fc) modulo 1, (baker map) (19) 

Mc(C,0 = + + + Q modulo 1, (cat map) (20) 



and 



^ ^ ^ 2^ sm (27r£), modulo 1, (standard map) (21) 



with < C, £ < 1. fee 2N is a parameter controlling the chaoticity of the map, i.e. the magnitude of the Lyapunov 
exponents. 

The resulting spectra are shown in Fig. [I] where we have also plotted the Lyapunov spectrum for elastic reflection 
as reference. To emphasize the conjugate pairs, the Lyapunov exponents are ordered as {A2jv-i+i> ^-2N+i}, with 
i = 1, . . . , 2N. Errors are estimated as in |20|| from the convergence of the exponents as a function of simulation time 
such that the time-dependent exponents did not deviate more than ±AA from their mean values during the second 
half of the simulation run. For high accuracy more than 10 7 disk-disk collisions and more than 5 • 10 6 disk-wall were 
simulated yielding errors less than ±0.001 for the exponents and less than ±0.002 for the pair sums. In the case of 
elastic reflection three Lyapunov exponents vanish. One exponent vanishes due to the neutral expansion behavior in 
the direction of the flow, a second due to the conservation of kinetic energy. The third exponents is zero due to the 
translational invariance of the system in the ^-direction pc| . The fourth vanishing exponent then verifies the pairing 
rule. For a hard disk system thermostated by deterministic scattering only two Lyapunov exponents vanish. The 
kinetic energy is now allowed to fluctuate around a mean value, so only the neutral expansion and the translational 
invariance remain. As we expect the maximum Lyapunov exponent increases with increasing chaoticity of the map, 
i.e. when going from a baker map (k = 2) to a cat map (k = 2) to a standard map (k = 100) (Results not plotted 
here show a similar behavior when k is increased for a given map.). The pairing rule for these models is satisfied with 
an error of ±0.002 in the pair sums. At this point we add a remark which might seem at first purely technical. For 
all simulations we used a symmetrical configuration, i.e. Eq. ([!]) is used for the upper wall whereas JVl and M.~ x are 
interchanged in Eq. ([!]) for the lower wall. Using the same scattering rules for both walls results in an asymmetry 
and eventually in an asymmetric Lyapunov spectrum even in equilibrium violating the pairing rule. 

B. NSS 

We move on to the nonequilibrium stationary state and turn first to the case of an imposed temperature gradient 
by the walls. Since in the thermodynamic limit the Lyapunov spectrum is mainly determined by the bulk behavior 
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we use in the following only a cat map with k = 2 as chaotic map M.. In order to determine the macroscopic state 
of the system the velocity and density profiles of the bulk are measured as well as the temperatures and velocities 
of the walls. Wall velocities are defined as the mean tangential velocity of the incoming and outgoing particles. 
Wall temperatures are defined as mean temperature of the incoming and outgoing fluxes, T = (Tj +T a ) /2, with 

Tj/ = (((v x — (v x ) x ) 2 ) x + [vy] y I [ u ! 7 1 ] y ) /2 where () x and [] y represent an average over the density p(v x ) and the 

flux <!> to and from the wall, respectively (see also Jjjf ). 



1. Heat flow 



Again, we first investigate a small system with four particles at n = 0.2 with high accuracy. In order to impose 
a temperature difference on the system we only have to choose two different parametrical temperatures T Ul T d (see 
Eq.(||) and p8|). Note that this also affects the derivatives of the collision matrix (Eqs. ( Jl7| , |l8|)). Figure |2| shows 
the spectra for this system under a temperature gradient, the numbers denoting the parametrical temperatures T u / d 
of the upper and the lower wall. In NSS the sum of the Lyapunov exponents is negative and exactly equal to the 
phase space contraction rate. T u — = 1 — 3 results in J] A; = —1.029 and T u — T d = 1 — 5 in A; = —2.703. 
We find again two vanishing Lyapunov exponents but with increasing temperature difference all nonzero exponents 
also increase in magnitude, the negative ones certainly stronger to yield an overall negative sum. The pair sums are 
also shown and, as we expect, the driving shifts the sums towards negative values thus destroying the symmetry. 
The deviations from the pairing rule are particularly strong for pairs with large i. Figure || shows the results for a 
36-particle system under the same setting. At least 2 • 10 6 disk-disk collisions and 2 • 10 5 disk-wall collisions have 
been simulated in each run. The spectra are plotted as connected lines only for graphical reasons, it is understood 
that the exponents are defined for integer i only. The change in the Lyapunov spectrum under thermal driving are 
similar to the four-particle system. Increasing the density from n — 0.2 to n — 0.6 results in a larger magnitude of 
all nonzero exponents due to the higher collision rate. The Kaplan- Yorke dimensionD^y [Fig. ^] is decreasing for 
increasing temperature gradient, with a larger dimensionality loss ADky for higher densities than for lower densities 
at given AT. As can immediately be guessed from the positive branch of the spectra thermal driving also results in 
an increasing Kolmogorov-Sinai entropy has [Fig. defined as the sum over all positive exponents, 

h K s= ]T Xi, (22) 

{A,>0} 

with increasing temperature gradient. So, as we expect from thermodynamics, thermal driving reduces the ordering 
of the system. Higher collision rates at higher densities lead to more viscous heating in the bulk and eventually to an 
increasing disorder (Hks) of the system. The second, lower data point at AT = shows Hks/N for elastic reflection 
as reference. 



2. Shear flow 



One way to model moving walls is to add some tangential momentum d to p x before and after the collision of a 
particle with the boundary (model I in [fl8||), 

(v f v f) _(S d oC+oS d ( P 1,pI), Pi > -d 

^Py)-\ Sd oC-oS d (p x ,pl), Px <-d, [Z6) 

with 



Sd(p i x ,pi) = (pi + d,p i y ). (24) 

In order to impose shear the shift d has only to be chosen with different signs for the upper and the lower wall. Note 
that these scattering rules are time reversible. Certainly, the drift also affects the derivatives of the collision matrix, 
Eqs. (|l7],|l8|), where p[ goes to p[ — d and p x to p x + d. Figure || shows the full Lyapunov spectra and the pair sums 
for a 36-particle system at n — 0.6 under shear while keeping T u — T d — 1 fixed. The negative exponents increase 
in magnitude whereas the positive branch changes very little. Before we take a closer look at the Kaplan- Yorke 
dimension and the Kolmogorov-Sinai entropy let us investigate another scattering rule (model III in Jl8|| ) which also 
models moving walls: 
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(p' x7 p' y ) = T- 1 oMoT4p x ,Pv) (25) 

with 

erf 

TAv^Vy) = I — ft ,exp(-p*/2T) I (26) 



Model III is still deterministic but no longer time reversible and only using this shear model a (numerical) equality 
between phase space contraction rate and entropy production was found in [fl8| . Note that p{. changes now to p{. — d 
and pi to p x — d in Eqs. ( |l7| , |l8| ). Figure || shows the corresponding Lyapunov spectra for a 36-particle system at 
ri = 0.6 under shear, but in contrast to model I the Lyapunov exponents of both the positive and the negative branch 
now increase in magnitude with increasing shear rate 7. For comparison the Lyapunov spectrum for the Chernov- 
Lebowitz model is also plotted, with Ek in /N and 7 equal to values obtained with model III at d = 1.5. Dky f° r 
both models is compared in Fig. [7(a)| and we see that the dimensionality loss with i ncrea sing shear rate 7 is stronger 
for model I than for model III. The graph of Hks-, plotted for both models in Fig. |7(b)| , asks for more explanation. 
Firstly, the overall behavior of hxs is increasing with larger 7, which seems to be the opposite of the observation 
made for the Chernov-Lebowitz shear model in j2Cj . But there the total kinetic energy is kept constant for all shear 
rates whereas here both models try to fix the wall temperature. Increasing shear results in an i ncre asing viscous heat 
production in the bulk which is reflected by a larger mean kinetic energy per particle [Fig. |3(a)[ |. Hence, the loss 
in hxs due to the ordering introduced by the shear is more than compensated by an increase of disorder due to a 
higher temperature in the bulk. Secondly, the only minor changes in the positive branch of the Lyapunov spectra for 
model I under shear yield an initially almost constant or even decreasing Kolmogorov-Sinai entropy. This, and the 



even more puzzling behavior of the wall temperature [Fig. 8(b)] can be explained by the fact that model I does not 
produce a Gaussian outgoing flux after the scattering. We found in |l8| that model I leads to an outgoing distribution 
with strong discontinuities in NSS whereas model III yields proper outgoing Gaussians. For comparison we have also 
computed Dky and hxs for the Chernov-Lebowitz model when it is approximately in the same macroscopic state as 
model III, i.e. we set the Chernov-Lebowitz system on the same kinetic energy shell and tried to find the appropriate 
shear parameter which re sults in the same shear rate. The Kaplan- Yorke dimension seems to be almost identical 



with that of model I, Fig. 7(a), and furthermore, the Kolmogorov-Sinai entropy now also increases with larger shear 



rate, only differing by a constant with the one obtained from model III, Fig. 7(b). This offset, depending on the 
special type of chaotic map chosen, originates from the fluctuating character of the model and should vanish in the 
thermodynamic limit. 



V. CONCLUSION 



We have calculated the full Lyapunov spectrum for a hard-disk fluid in equilibrium and nonequilibrium steady states 
thermostated by deterministic scattering. Since the model allows for fluctuations around a mean total energy only 
two vanishing Lyapunov exponents are found in both equilibrium and nonequilibrium states. In nonequilibrium the 
system is dissipative with a mean phase-space contraction rate smaller than zero. The magnitude of the Lyapunov 
exponents increases with increasing temperature gradient or shear rate, with a stronger increase for the negative 
branch. Thus both heat and shear flow situations result in a decreasing Kaplan- Yorke dimension and an increasing 
Kolmogorov-Sinai entropy with stronger nonequilibrium. Due to the inhomogeneous driving at the boundary the 
pairing rule does not hold. We did not verify the relation between Lyapunov exponents and transport coefficients (as 
e.g. in p2| ) since in absence of a pairing rule this would be equivalent to checking the relation between phase space 
contraction and entropy production rates, where the latter has been done in fjsj] . The main difference of our shear 
flow model III and the Chernov-Lebowitz shear flow model is the fact that in the latter the total kinetic energy is 
fixed whereas our models tries to fix the wall temperature. A direct comparison of both models reveals that they yield 
identical Kaplan- Yorke dimensions and only differ by a constant in the Kolmogorov-Sinai entropy, which vanishes in 
the thermodynamic limit. 
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I 

FIG. 1. Lyapunov spectra for a four particle system with density n — 0.2 in equilibrium. The solid line shows the spectrum 
for elastic reflection, i.e. the identity map is used instead of a chaotic one. The other spectra are obtained by using a baker 
map with k — 2, a cat map with k — 2 and a standard map with k — 100, see Eqs. 



it her spe c 




FIG. 2. Lyapunov spectra for a four particle system at density n — 0.2 in NSS. The imposed temperatures for the upper and 
the lower wall are indicated by the numbers. The respective pair sums (scaled by a factor 1/2 for graphical reasons) are also 
plotted near the middle line. 
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FIG. 3. Full Lyapunov spectra for a 36-particle system under an imposed temperature gradient at number densities n — 0.2 
and n = 0.6. The numbers indicate the parametrical temperatures T u , Td- The respective pair sums (scaled by a factor 1/2 
for graphical reasons) are also plotted near the middle line. 
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FIG. 4. The Kaplan- Yorke dimension Dry and the Kolmogorov-Sinai entropy has per particle for a 36-particle system 
under a temperature gradient at number densities n = 0.2 and n = 0.6. AT denotes the measured temperature difference 
between the upper and the lower wall. The second, lower data point at AT = gives hxs/N for elastic reflection as reference. 
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FIG. 5. Full Lyapunov spectra for a 36-particle system under shear at number density n = 0.6, model I 
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FIG. 6. Full Lyapunov spectra for a 36-particle system under shear at number density n = 0.6, model III. The spectrum for 
the Chernov-Lebowitz model results from a simulation with the same mean kinetic energy per particle and approximately the 
same shear rate as model III with T u — T d — 1.0 and d = 1.5. 
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FIG. 7. The Kaplan- Yorke dimension Dry and the Kolmogorov- Sinai entropy has per particle for a 36-particle system 
under shear at number density n — 0.6 as a function of the shear rate 7 (model I, model III and Chernov- Lebowitz model). 
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FIG. 8. Kinetic energy per particle and measured wall temperature for a 36-particle system under shear at number density 
n = 0.6 as a function of the shear rate 7 (model I and model III). 
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